COVID-19 Data Analysis

Rt per provincia

Rt di COVID-19 stimato per le province italiane.

Max Pierini, Alessio Pamovio & NOTIZIÆ Telegram Channel


NB: gli $R_t$ stimati da EpiDataItalia sono parametrizzati attualmente sui soli nuovi casi riportati al Dipartimento di Protezione Civile e non sono sovrapponibili al numero di riproduzione stimato da ISS-EpiCentro che utilizza un modello EpiEstim (A new framework and software to estimate time-varying reproduction numbers during epidemics, Cori-Ferguson-Fraser 2013) per il quale sono necessari dati attualmente non disponibili al pubblico (il numero di casi sintomatici con data di inizio sintomi, il numero di casi sintomatici importati da un’altra regione o dall’estero). Per info sul modello in uso da ISS-EpiCentro consultare la pagina FAQ sul calcolo del Rt di ISS. Per chiedere al Governo l'accesso pubblico a tutti i dati grezzi disaggregati di COVID-19 in Italia consigliamo l'adesione alla Petizione DatiBeneComune.

Method

A simple method is presented to estimate effective reproduction number $R_t$ of COVID-19 in italian provinces with a Markov chain Monte Carlo and Poisson likelihood parametrized on daily new cases.

Details

New cases $y_r$ for each $r$ italian province (source Dipartimento Protezione Civile), will be smoothed with rolling mean (gaussian, window 7, std 2). Smoothed new cases will be adjusted to be $>0$ to avoid negative values (due to data error corrections).

For each day $t$, smoothed new cases $y_{r,t}$ will be supposed distributed as Poisson with $\lambda_{r,t}$ parameter

$$ y_{r,t} \sim \mathcal{P}(\lambda_{r,t}) $$

where $\lambda_{r,t}$ is defined by the serial interval inverse $\gamma$, previous day smoothed new cases $k_{r,t-1}$ and effective reproduction number in time $R_{r,t}$ (ref: Bettencourt & Ribeiro 2008)

$$ \lambda_{r,t} = k_{r,t-1} e^{\gamma (R_{r,t} - 1)} $$

The serial interval SI is supposed to be distributed as Gamma, with mean $\mu=7.5$ and standard deviation $\sigma=3.4$ (ref: Li, Ghuan et Al. 2020a)

$$ \mathbf{SI} \sim \Gamma(\mu_{=7.5}, \sigma_{=3.4}) $$

so that $\gamma$ is distributed as Inverse Gamma

$$ \gamma \sim \Gamma^{-1}(\mu_{=7.5}, \sigma_{=3.4}) $$

Parameters $R_{{p,t}}$, for each province $p$ and time $t$ (day) will be distributed as half normal with mean equal to previous day posteriors $R_{{p,t-1}}$ and precision $\tau$

$$ R_{{p,t}} \sim \mathcal{{N}}^+(R_{{p,t-1}}, \tau) $$

where, first day $R_{{p,0}}$ (outcome) is set to zero

$$ R_{{p,0}} = 0 $$

and $\tau$ is defined as $1 / \sigma^2$, where $\sigma$ is the standard deviation of $R_t$ priors, previuosly estimated for italian provinces (using the same method but with overarching uninformative prior distribution of $\tau$) as

$$ \sigma \sim \mathcal{N}(\mu_{=0.299}, \sigma_{=0.004}) $$

If previous new cases are zero $k_{p,t-1}=0$, parameter $R_{p,t}$ is undefined, given the chosen function for $\lambda_{p,t}$ parameter of Poisson likelihood, even if it should be $R_{p,t}=0$ (no new cases means null effective reproduction number). Thus, in these cases, priors of $R_{p,t}$ will be forced to

$$ R_{p,t} \sim \mathcal{N}^+(0, \tau) \;,\; k_{p,t-1}=0 $$

and previous new cases will be forced to $k_{p,t-1}=1$, so that $\lambda_{p,t}$ will be

$$ \lambda_{p,t} = e^{\gamma( \mathcal{N}^+(0, \tau) - 1 )} \;,\; k_{p,t-1}=0 $$

A Markov chain Monte Carlo will be used with Metropolis-Hasting algorithm and Gibbs sampler (adapt 500, warmup 1000, samples 1000, chains 4, thin 1) in Python 3.8.5 with pyjags==1.3.7 and JAGS 4.3.0.

model {
    # Overarching Rt standard deviation
    sigma_R ~ dnorm( 0.29929894605552754 , 1 / 0.004218809765237365 )
    tau_R <- 1 / sigma_R^2


    # Serial interval distribution
    SI_mu <- 7.5
    SI_sd <- 3.4
    SI_sh <- SI_mu^2 / SI_sd^2
    SI_ra <- SI_mu / SI_sd^2
    SI ~ dgamma( SI_sh , SI_ra )
    gamma <- 1 / SI

    # First Rt prior
    R[1] <- 0
    for ( t in 2:T ) {
        # Rt prior for k>0
        Rpr[t] ~ dnorm( R[t-1] , tau_R )  T(0,)
        # Rt prior for k=0
        Rnu[t] ~ dnorm( 0 , tau_R )  T(0,)

        # Define Rt prior
        R[t] <- ifelse( k[t-1]==0 , Rnu[t] , Rpr[t] )
        # Avoid k=0 (undefined Rt)
        K[t] <- ifelse( k[t-1]==0, 1 , k[t-1] )

        # Poisson likelihood
        lambda[t] <- K[t] * exp( gamma * ( R[t] - 1 ) )
        y[t] ~ dpois( lambda[t] )
    }
}
DONE!

Time series

Latest Rt (IQR)

Latest Rt (HPDI)

Latest Rt (MAP)


-> 2020-11-10 18:34:32.988763
<- 2020-11-10 19:53:20.031080
Completed in 1:18:47.042317

© 2020 Max Pierini. Thanks to Sandra Mazzoli & Alessio Pamovio

Exported from Italia/Rt_Province.ipynb committed by maxdevblock on Tue Nov 10 20:00:32 2020 revision 4, 9600503